set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)

  data <- read.dta('Netherlands (2012).dta')

	model <- clmm(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + Perc_seats + Prime_minister * Perc_seats + Cabinet * Perc_seats + Opposition*Perc_seats + (1 | Party) + (1 | pid), data = data, HESS = TRUE, link = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
	vce <- vcov(model)
	vce <- vce[-c(12:13), ]      
	vce <- vce[, -c(12:13)]       
	
	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

	partHi <- c()
	partLo <- c()
	partMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()


	for(i in 0:50){
	## pm party 

		agg <- simCoef[, 5] * 1 + 
				simCoef[, 6] * 0 + 
				simCoef[, 7] * 0 + 
				simCoef[, 8] * i + 
				simCoef[, 9] * i + 
				simCoef[, 10] * 0 + 
				simCoef[, 11] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		pmHi[i + 1] <- quantile(pred, 0.975)
		pmLo[i + 1] <- quantile(pred, 0.025)
		pmMed[i + 1] <- quantile(pred, 0.5)


	## partner party 
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 1 + 
	  simCoef[, 7] * 0 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * i + 
	  simCoef[, 11] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partHi[i + 1] <- quantile(pred, 0.975)
		partLo[i + 1] <- quantile(pred, 0.025)
		partMed[i + 1] <- quantile(pred, 0.5)

	## opp party 
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 1 + 
	  simCoef[, 8] * i + 
	  simCoef[, 9] * 0 + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
		oppLo[i + 1] <- quantile(pred, 0.025)
		oppMed[i + 1] <- quantile(pred, 0.5)

  ## No seats
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 0 + 
    simCoef[, 7] * 0 + 
    simCoef[, 8] * i + 
    simCoef[, 9] * 0 + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * 0 
      
      pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
      pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
      pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
      pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
      pred5 <- 1 - pred1 - pred2 - pred3 - pred4
      
      pred <- pred5 * 5 + 
        pred4 * 4 + 
        pred3 * 3 + 
        pred2 * 2 + 
        pred1
      
      otherHi[i + 1] <- quantile(pred, 0.975)
      otherLo[i + 1] <- quantile(pred, 0.025)
      otherMed[i + 1] <- quantile(pred, 0.5)
    
    }


## now paint a pretty picture.... ##

	seats <- seq(0, 50, 1)        
  	
	#Now create the figure for Denmark
	par(mar=c(3.6,4,1,2))      
	
	plot(1, 1, type = 'n',
		xlim = c(0, 50),
		ylim = c(1, 5),
		xlab = '',
		ylab = '',
		main = '',
		axes = FALSE)
	axis(1,at=seq(0,50,25),labels=T)          
	axis(2)
	lines(seats, pmMed, col = rgb(1, 0, 0, 0.2), lwd = 3)

## darken bands in high density range ##
	seatHi <- quantile(data$Perc_seats[data$Prime_minister == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$Perc_seats[data$Prime_minister == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], pmMed[seats >= seatLo & seats <= seatHi], col = rgb(1, 0, 0, 0.3), lwd = 15)
	
	lines(seats, partMed, col = rgb(0, 0, 1, 0.2), lwd = 3)

	seatHi <- quantile(data$Perc_seats[data$Cabinet == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$Perc_seats[data$Cabinet == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], partMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 1, 0.3), lwd = 15)

	lines(seats, oppMed, col = rgb(0, 0, 0, 0.2), lwd = 3)

	seatHi <- quantile(data$Perc_seats[data$Opposition == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$Perc_seats[data$Opposition == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], oppMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 0, 0.3), lwd = 15)

  points(0, y = otherMed[1], col = rgb(1, 0, 1, 0.3), lwd = 15)  

	mtext("Responsibility attribution", side=2, line=2.5, cex=1.5)        
	mtext("Perceived seat %", side=1, line=2.2, cex=1.5)                  


############
###LEGEND###
############

legend('bottomright', legend = c('PM', 'Partner', 'Opposition', 'No seats'), fill = c(rgb(1, 0, 0, 0.3),	
	rgb(0, 0, 1, 0.3), rgb(0, 0, 0, 0.3), rgb(1, 0, 1, 0.3)), bty = 'n', border = NA, , cex=1.25)

###################
###EXPORT FIGURE###
###################

dev.copy(png,'Figure 1b NL 2012 - party types over seat share.png', height = 490, width = 633)	
dev.off()  
